Lendo o Shapefile crime_mg:
crime_mg<- readOGR(dsn = "crime_mg", layer = "crime_mg",verbose = TRUE, use_iconv = TRUE, p4s = "+proj=longlat +ellps=WGS84")
OGR data source with driver: ESRI Shapefile
Source: "/Users/ricardo/Tresors/zz-pessoal/FGV/git/Trabalhos/GAEE/Tarefa 4/crime_mg", layer: "crime_mg"
with 754 features
It has 17 fields
Integer64 fields read as strings: POP_RUR POP_URB POP_FEM POP_MAS
names(crime_mg)
[1] "CODMUNI" "ID" "MUNIC" "AREA" "INDICE94" "INDICE95" "GINI_91" "POP_94" "POP_RUR" "POP_URB"
[11] "POP_FEM" "POP_MAS" "POP_TOT" "URBLEVEL" "PIB_PC" "X_COORD" "Y_COORD"
summary(crime_mg)
Object of class SpatialPolygonsDataFrame
Coordinates:
min max
x -51.06258 -39.85724
y -22.91696 -14.23725
Is projected: FALSE
proj4string : [+proj=longlat +ellps=WGS84]
Data attributes:
CODMUNI ID MUNIC AREA INDICE94 INDICE95
Min. : 10 Min. : 0.0 Abadia dos Dourados: 1 Min. : 40.3 Min. : 0.260 Min. : 0.420
1st Qu.:1852 1st Qu.:188.2 Abaeté : 1 1st Qu.: 208.9 1st Qu.: 8.143 1st Qu.: 9.643
Median :3645 Median :377.5 Abre-Campo : 1 Median : 371.4 Median :12.060 Median :13.975
Mean :3636 Mean :377.5 Acaiaca : 1 Mean : 779.4 Mean :13.329 Mean :15.449
3rd Qu.:5444 3rd Qu.:566.8 Açucena : 1 3rd Qu.: 827.5 3rd Qu.:17.203 3rd Qu.:19.820
Max. :7220 Max. :755.0 Água Boa : 1 Max. :13292.1 Max. :41.300 Max. :47.690
(Other) :748
GINI_91 POP_94 POP_RUR POP_URB POP_FEM POP_MAS POP_TOT
Min. :0.0000 Min. : 820 0 : 33 0 : 33 0 : 33 0 : 33 Min. : 0
1st Qu.:0.5129 1st Qu.: 4724 1219 : 3 1374 : 3 1042 : 2 1226 : 2 1st Qu.: 4295
Median :0.5578 Median : 8602 1994 : 3 10429 : 2 1070 : 2 1483 : 2 Median : 8216
Mean :0.5330 Mean : 21640 4995 : 3 1120 : 2 1368 : 2 1539 : 2 Mean : 20865
3rd Qu.:0.5960 3rd Qu.: 18054 12472 : 2 11996 : 2 1449 : 2 1643 : 2 3rd Qu.: 17710
Max. :0.7127 Max. :2079280 1252 : 2 1257 : 2 1619 : 2 1658 : 2 Max. :2020161
(Other):708 (Other):710 (Other):711 (Other):711
URBLEVEL PIB_PC X_COORD Y_COORD
Min. :0.0000 Min. : 0 Min. :-50.81 Min. :-22.81
1st Qu.:0.3743 1st Qu.: 1665 1st Qu.:-45.55 1st Qu.:-21.19
Median :0.5445 Median : 2446 Median :-44.06 Median :-20.01
Mean :0.5373 Mean : 3036 Mean :-44.22 Mean :-19.81
3rd Qu.:0.7120 3rd Qu.: 3525 3rd Qu.:-42.76 3rd Qu.:-18.77
Max. :0.9970 Max. :37728 Max. :-40.03 Max. :-14.46
Mapa de Minas Gerais com os municípios, como no shapefile, sem tema:
tmap::qtm(crime_mg,title = "Mapa de Minas Gerais")
Linking to GEOS 3.6.1, GDAL 2.1.3, proj.4 4.9.3
Calculo de Moran’s I para verificação da auto-correlação espacial das variáveis. Aqui usamos a metodologia Rainha (queen) na matriz de vizinhança:
crime_mg_nb = poly2nb(crime_mg, queen=TRUE, row.names=crime_mg$X_COORD)
crime_mg_w <- nb2listw(crime_mg_nb, style="W")
cmg_munic <- as.numeric(crime_mg$MUNIC)
cmg_area <- as.numeric(crime_mg$AREA)
cmg_indice94 <- as.numeric(crime_mg$INDICE94)
cmg_indice95 <- as.numeric(crime_mg$INDICE95)
cmg_gini_91 <- as.numeric(crime_mg$GINI_91)
cmg_pop_94 <- as.numeric(crime_mg$POP_94)
cmg_pop_rur <- as.numeric(crime_mg$POP_RUR)
cmg_pop_urb <- as.numeric(crime_mg$POP_URB)
cmg_pop_fem <- as.numeric(crime_mg$POP_FEM)
cmg_pop_mas <- as.numeric(crime_mg$POP_MAS)
cmg_pop_tot <- as.numeric(crime_mg$POP_TOT)
cmg_urblevel <- as.numeric(crime_mg$URBLEVEL)
cmg_pib_pc <- as.numeric(crime_mg$PIB_PC)
moran_i_munic <- moran(cmg_munic,crime_mg_w, length(crime_mg_nb), Szero(crime_mg_w))
moran_i_area <- moran(cmg_area,crime_mg_w, length(crime_mg_nb), Szero(crime_mg_w))
moran_i_indice94 <- moran(cmg_indice94,crime_mg_w, length(crime_mg_nb), Szero(crime_mg_w))
moran_i_indice95 <- moran(cmg_indice95,crime_mg_w, length(crime_mg_nb), Szero(crime_mg_w))
moran_i_gini_91 <- moran(cmg_gini_91,crime_mg_w, length(crime_mg_nb), Szero(crime_mg_w))
moran_i_pop_94 <- moran(cmg_pop_94,crime_mg_w, length(crime_mg_nb), Szero(crime_mg_w))
moran_i_pop_rur <- moran(cmg_pop_rur,crime_mg_w, length(crime_mg_nb), Szero(crime_mg_w))
moran_i_pop_urb <- moran(cmg_pop_urb,crime_mg_w, length(crime_mg_nb), Szero(crime_mg_w))
moran_i_pop_fem <- moran(cmg_pop_fem,crime_mg_w, length(crime_mg_nb), Szero(crime_mg_w))
moran_i_pop_mas <- moran(cmg_pop_mas,crime_mg_w, length(crime_mg_nb), Szero(crime_mg_w))
moran_i_pop_tot <- moran(cmg_pop_tot,crime_mg_w, length(crime_mg_nb), Szero(crime_mg_w))
moran_i_urblevel <- moran(cmg_urblevel,crime_mg_w, length(crime_mg_nb), Szero(crime_mg_w))
moran_i_pib_pc <- moran(cmg_pib_pc,crime_mg_w, length(crime_mg_nb), Szero(crime_mg_w))
Mostrando todas as auto-correlações:
moran <- c("label","i")
label <- c("moran_i_munic",
"moran_i_area",
"moran_i_indice94",
"moran_i_indice95",
"moran_i_gini_91",
"moran_i_pop_94",
"moran_i_pop_rur",
"moran_i_pop_urb",
"moran_i_pop_fem",
"moran_i_pop_mas",
"moran_i_pop_tot",
"moran_i_urblevel",
"moran_i_pib_pc"
)
moran_i <- c(moran_i_munic$I
,moran_i_area$I
,moran_i_indice94$I
,moran_i_indice95$I
,moran_i_gini_91$I
,moran_i_pop_94$I
,moran_i_pop_rur$I
,moran_i_pop_urb$I
,moran_i_pop_fem$I
,moran_i_pop_mas$I
,moran_i_pop_tot$I
,moran_i_urblevel$I
,moran_i_pib_pc$I
)
moran <- data.frame(label = label, moran_i = moran_i)
moran[order(moran$moran_i,decreasing = TRUE),]
Mostrando Moran’s I das variáveis AREA, INDICE94 e INDICE95, que possuem os maiores I, e POP_URB, que possui o menor:
{
moran.plot(x = cmg_area,listw = crime_mg_w,labels = FALSE)
title("Moran's I de AREA")
moran.plot(x = cmg_indice94,listw = crime_mg_w,labels = FALSE)
title("Moran's I de INDICE94")
moran.plot(x = cmg_indice95,listw = crime_mg_w,labels = FALSE)
title("Moran's I de INDICE95")
moran.plot(x = cmg_pop_urb,listw = crime_mg_w,labels = FALSE)
title("Moran's I de POP_URB")
}
Calculando LISA
Verificando a média de links entre vizinhos:
crime_mg_nb
Neighbour list object:
Number of regions: 754
Number of nonzero links: 4302
Percentage nonzero weights: 0.7567069
Average number of links: 5.70557
Como a média de links é 5.7, passamos este valor como parâmetro com o comando “mean(card(crime_mg_nb))” para o cálculo do LISA para as variáveis AREA - INDICE94 - INDICE95 - POP_URB:
crime_mg$X_COORD
[1] -47.49321 -45.38105 -42.43796 -43.10176 -42.40581 -42.27959 -48.09214 -45.44244 -40.96760 -41.63569
[11] -41.19622 -44.66089 -44.68020 -46.62777 -42.75312 -46.01826 -43.72144 -40.71614 -42.00387 -46.32413
[21] -46.19588 -41.93930 -43.41420 -41.67652 -43.14477 -43.38958 -42.80526 -46.58660 -44.29237 -43.76841
[31] -42.88019 -42.16082 -44.21976 -43.41933 -41.93442 -48.26069 -44.23796 -42.49868 -49.11842 -46.12565
[41] -45.19000 -46.99052 -46.95313 -45.52037 -46.14603 -42.84672 -46.00288 -42.88961 -41.26131 -44.23362
[51] -44.84372 -43.86438 -45.99177 -40.57566 -46.39785 -43.52186 -42.28686 -43.80915 -43.07197 -43.95786
[61] -43.12605 -43.48852 -43.96807 -42.43910 -44.04636 -42.57421 -40.60139 -44.20018 -43.78513 -43.11836
[71] -45.51844 -45.66553 -44.50138 -43.58578 -45.30816 -44.11943 -46.54138 -43.47979 -42.37122 -46.19729
[81] -44.80967 -44.22466 -46.24602 -46.17535 -46.42886 -42.96351 -43.23695 -44.58306 -45.62654 -42.70726
[91] -44.16087 -46.36739 -44.09327 -46.59474 -45.17579 -46.39264 -44.47074 -45.80338 -41.49338 -49.47046
[101] -44.39202 -43.63525 -41.92847 -42.76596 -46.35720 -45.15727 -46.08641 -46.06212 -45.26101 -41.72284
[111] -45.42651 -46.24080 -49.75845 -45.25514 -45.83823 -48.65129 -46.19711 -45.77307 -45.19620 -42.63518
[121] -49.24278 -45.29543 -41.90053 -43.62184 -42.47477 -47.03152 -44.16349 -49.56424 -41.83155 -43.66496
[131] -46.14776 -42.26532 -41.59302 -43.73197 -43.83725 -42.11535 -42.16643 -43.05118 -45.67330 -40.85270
[141] -43.16063 -45.21992 -44.89192 -45.17202 -44.71230 -46.18546 -46.12962 -44.66230 -50.80962 -44.62952
[151] -45.84430 -44.49512 -43.95010 -47.87196 -46.94624 -42.66952 -43.50941 -41.49439 -44.96505 -45.68552
[161] -41.31457 -49.16922 -43.23795 -41.66556 -42.39806 -43.02568 -43.38604 -47.24249 -44.22839 -44.76685
[171] -42.80558 -42.83844 -49.09380 -41.78283 -46.22109 -44.49410 -48.35280 -45.41821 -41.72209 -43.50147
[181] -44.88534 -45.10409 -45.79407 -46.06016 -43.85135 -43.71929 -47.63420 -43.80803 -41.38163 -45.92734
[191] -44.08560 -45.42834 -44.41458 -44.21446 -45.66895 -44.59043 -42.26304 -47.13203 -42.72135 -42.20463
[201] -43.29637 -44.20972 -45.96833 -46.00674 -42.45298 -43.41827 -45.56686 -42.80937 -43.82191 -45.27968
[211] -44.36035 -46.68877 -44.79921 -44.43959 -43.64483 -45.29302 -46.79451 -42.97382 -44.27892 -43.51764
[221] -43.65886 -43.17408 -42.68796 -43.00446 -42.16505 -41.49780 -42.61690 -44.93895 -46.25656 -40.93218
[231] -42.09658 -43.25840 -42.91237 -45.15475 -42.81011 -43.99833 -42.93293 -45.54055 -43.16655 -45.89812
[241] -47.62720 -41.82026 -45.60380 -42.02492 -43.99986 -42.25049 -44.10489 -42.60332 -44.32296 -41.90697
[251] -43.01434 -45.99563 -46.02793 -42.49148 -45.79619 -47.75566 -42.25461 -43.57525 -46.28983 -45.82582
[261] -42.04826 -43.23039 -40.74375 -44.94769 -42.11291 -42.95870 -42.35518 -44.43566 -45.54879 -46.10599
[271] -46.76576 -44.53501 -42.26058 -44.28552 -43.45775 -41.51931 -41.84971 -49.16301 -40.85224 -49.01385
[281] -44.07496 -41.49296 -45.83345 -42.49620 -43.88004 -41.93332 -42.95252 -47.77042 -42.81914 -45.92636
[291] -43.03085 -46.80881 -43.05771 -43.01696 -47.12282 -46.69458 -42.78351 -46.74889 -42.69890 -49.85956
[301] -45.54329 -42.28752 -43.96064 -46.58465 -44.80533 -47.11019 -44.11867 -46.40874 -44.79231 -44.87152
[311] -44.29277 -44.72349 -45.76220 -44.93634 -45.80854 -46.29480 -47.86486 -44.95249 -41.96208 -44.42951
[321] -44.28073 -42.40428 -41.67196 -42.60780 -49.92467 -46.15867 -47.46733 -43.31345 -41.24347 -43.80298
[331] -43.34578 -44.19357 -44.53723 -41.66136 -45.40179 -42.83254 -42.84730 -41.86216 -43.36308 -47.04784
[341] -44.75876 -44.91725 -41.82427 -41.51701 -49.44074 -45.10546 -46.23148 -44.46077 -46.77563 -44.58506
[351] -43.61058 -41.67828 -41.11247 -49.52466 -44.83189 -50.37773 -44.71800 -43.72156 -40.31645 -46.75804
[361] -46.60800 -42.72727 -43.69300 -41.75105 -43.44411 -44.95351 -45.53340 -44.04658 -42.63187 -44.46296
[371] -43.99919 -41.04881 -41.03704 -42.72650 -43.16084 -45.97010 -44.12962 -40.31486 -44.34806 -43.43648
[381] -46.52044 -41.84221 -46.75294 -45.49371 -44.65812 -44.06599 -46.29684 -46.50148 -43.91002 -41.60307
[391] -45.38542 -43.49854 -42.45140 -44.70590 -45.03210 -45.03551 -42.68481 -44.32724 -43.88068 -50.67722
[401] -44.28486 -44.91229 -45.70727 -45.91916 -44.34937 -42.06419 -42.95220 -44.17779 -42.04783 -41.86315
[411] -41.10835 -43.02729 -44.67303 -45.31263 -43.34615 -42.07090 -42.96997 -42.64959 -45.16636 -45.20746
[421] -40.70877 -43.03168 -44.44209 -43.31573 -43.76842 -41.93533 -42.31006 -42.99897 -44.06631 -46.05403
[431] -40.73267 -46.37550 -41.50532 -41.39193 -43.34962 -42.59135 -42.44463 -44.59889 -44.14048 -42.39810
[441] -42.63601 -44.00836 -45.40095 -44.02541 -45.50091 -44.45995 -48.89517 -43.15817 -46.34118 -47.49709
[451] -46.96555 -46.57630 -43.94379 -42.45094 -45.40588 -44.64335 -43.40523 -46.31461 -42.43851 -41.47467
[461] -46.53731 -42.18540 -40.50734 -45.49600 -44.60592 -45.28030 -43.01226 -43.89931 -41.53249 -47.71389
[471] -46.43064 -44.99744 -43.59176 -41.97493 -43.98354 -45.28965 -44.69688 -43.53341 -44.74269 -43.69686
[481] -46.38547 -43.66598 -41.25998 -41.51534 -45.57246 -45.69416 -43.43626 -42.31288 -40.36231 -44.70464
[491] -44.61044 -46.87930 -45.76612 -45.82705 -44.45277 -44.95725 -44.48098 -44.27518 -43.19171 -46.63207
[501] -46.41648 -47.06932 -42.25492 -42.98432 -42.86879 -41.09830 -42.54057 -41.18228 -42.72473 -45.23163
[511] -42.16079 -45.46335 -44.36013 -47.52156 -44.06279 -43.74113 -43.13725 -44.64038 -45.05287 -47.15188
[521] -45.07785 -41.56500 -43.31730 -42.72183 -44.17341 -44.25210 -45.83891 -44.42691 -48.65832 -43.30649
[531] -45.50943 -45.59827 -42.37189 -44.85914 -43.03453 -44.87062 -46.06173 -48.64491 -46.01222 -46.54948
[541] -41.54432 -44.93298 -42.89923 -43.03143 -43.08218 -41.80911 -45.95200 -44.94014 -44.09348 -48.91177
[551] -46.85789 -46.40331 -43.15821 -44.08090 -43.59392 -46.37173 -44.12725 -45.60957 -43.88958 -43.78210
[561] -42.35612 -42.45291 -44.29963 -41.16264 -43.77735 -45.96197 -42.99298 -44.06237 -45.09436 -43.77328
[571] -42.66293 -40.55101 -42.90916 -43.50816 -44.35619 -43.16596 -46.31690 -42.41469 -43.17127 -43.18086
[581] -43.80368 -43.09106 -44.37965 -43.04526 -42.84922 -47.55183 -42.20835 -40.51328 -43.79208 -43.06769
[591] -47.23560 -42.18431 -40.02573 -43.46690 -42.10530 -43.53213 -42.82593 -42.39977 -45.55776 -47.51220
[601] -43.86069 -42.26607 -43.05533 -40.12327 -42.36111 -46.27435 -42.11675 -43.95672 -41.39844 -44.08352
[611] -45.68169 -46.01767 -50.33390 -45.50389 -42.55058 -43.92583 -43.19073 -44.06586 -45.07629 -41.89212
[621] -42.55613 -43.64855 -43.68060 -44.94793 -42.81468 -42.62218 -43.26908 -40.25531 -45.27807 -43.24890
[631] -44.16770 -43.50954 -45.07244 -43.97797 -42.90958 -45.07765 -44.99714 -49.88219 -42.28881 -42.82680
[641] -42.30198 -45.64853 -44.82851 -43.32174 -43.36258 -45.59533 -45.98589 -46.44325 -45.95571 -43.92105
[651] -44.29027 -42.16253 -41.13876 -42.19143 -41.85145 -42.79912 -43.02130 -43.97927 -42.14351 -44.56289
[661] -45.54713 -41.37432 -42.67596 -42.69985 -41.77816 -45.03896 -42.71632 -46.66230 -42.59288 -42.58146
[671] -45.42453 -46.50826 -45.77652 -42.55856 -45.03375 -47.02728 -43.24157 -45.02241 -44.57803 -47.14284
[681] -44.97449 -44.47504 -45.84660 -42.40123 -46.22298 -42.91781 -43.11113 -46.15979 -43.25835 -43.35847
[691] -43.10017 -43.62730 -42.46251 -44.48242 -43.22733 -45.77631 -46.69787 -40.27975 -46.10095 -44.55666
[701] -43.39589 -44.30680 -43.20498 -45.80846 -43.30811 -42.05228 -42.15477 -45.04837 -43.26382 -41.98139
[711] -46.89534 -46.16463 -43.69368 -41.91384 -42.86211 -41.37456 -42.61551 -44.16510 -45.84037 -43.03073
[721] -46.39741 -42.06352 -45.23606 -45.04473 -45.49229 -41.64328 -48.77868 -42.83191 -45.80486 -42.97899
[731] -44.89385 -42.05144 -47.97703 -48.36030 -40.66907 -46.82458 -42.73167 -45.46964 -46.33885 -45.40760
[741] -44.73364 -43.94599 -46.88179 -45.41862 -48.35473 -43.93210 -42.88370 -42.28546 -42.31654 -45.12150
[751] -42.66643 -42.31988 -42.85113 -42.55222
Plot dos mapas LISA:
crime_mg$LISA_AREA_p <- LISA_AREA$p
crime_mg$LISA_INDICE94_p <- LISA_INDICE94$p
crime_mg$LISA_INDICE95_p <- LISA_INDICE95$p
crime_mg$LISA_POP_URB_p <- LISA_POP_URB$p
crime_mg$LISA_AREA_corr <- LISA_AREA$correlation
crime_mg$LISA_INDICE94_corr <- LISA_INDICE94$correlation
crime_mg$LISA_INDICE95_corr <- LISA_INDICE95$correlation
crime_mg$LISA_POP_URB_corr <- LISA_POP_URB$correlation
tmap::tm_shape(crime_mg, simplify = 1) +
tmap::tm_polygons() +
tmap::tm_shape(crime_mg, simplify = 1) +
tmap::tm_fill(c("LISA_AREA_p","LISA_AREA_corr"), midpoint = 0) +
tmap::tm_style("natural")
tmap::tm_shape(crime_mg) +
tmap::tm_polygons() +
tmap::tm_shape(crime_mg, simplify = 1) +
tmap::tm_fill(c("LISA_INDICE94_p","LISA_INDICE94_corr"), midpoint = 0) +
tmap::tm_style("natural")
tmap::tm_shape(crime_mg) +
tmap::tm_polygons() +
tmap::tm_shape(crime_mg, simplify = 1) +
tmap::tm_fill(c("LISA_INDICE95_p","LISA_INDICE95_corr"), midpoint = 0) +
tmap::tm_style("natural")
tmap::tm_shape(crime_mg) +
tmap::tm_polygons() +
tmap::tm_shape(crime_mg, simplify = 1) +
tmap::tm_fill(c("LISA_POP_URB_p","LISA_POP_URB_corr"), midpoint = 0) +
tmap::tm_style("natural")
Fizemos também a implementação dos mapas LISA no GeoDa e estes foram diferentes do apresentado aqui. Abaixo, os gráficos gerados no GeoDa:
knitr::include_graphics("MoransI-Geoda-Area.png")
knitr::include_graphics("LISA Significance - Geoda - Area.png")
knitr::include_graphics("LISA Cluster - Geoda - Area.png")
knitr::include_graphics("MoransI-Geoda-INDICE94.png")
knitr::include_graphics("LISA Significance - Geoda - INDICE94.png")
knitr::include_graphics("LISA Cluster - Geoda - INDICE94.png")
knitr::include_graphics("MoransI-Geoda-INDICE95.png")
knitr::include_graphics("LISA Significance - Geoda - INDICE95.png")
knitr::include_graphics("LISA Cluster - Geoda - INDICE95.png")
knitr::include_graphics("MoransI-Geoda-POP_URB.png")
knitr::include_graphics("LISA Significance - Geoda - POP_URB.png")
knitr::include_graphics("LISA Cluster - Geoda - POP_URB.png")
Regressão linear simples:
ap <- as.data.frame(cbind(crime_mg$MUNIC, crime_mg$AREA, crime_mg$INDICE94,
crime_mg$INDICE95, crime_mg$GINI_91, crime_mg$POP_94,
crime_mg$POP_RUR, crime_mg$POP_URB, crime_mg$POP_FEM,
crime_mg$POP_MAS, crime_mg$POP_TOT, crime_mg$URBLEVEL,
crime_mg$PIB_PC))
colnames(ap) <- c("ID",names(crime_mg@data[4:15]))
head(ap)
lmK <- lm(formula = crime_mg$INDICE95 ~ URBLEVEL, data = ap)
summary(lmK)
Call:
lm(formula = crime_mg$INDICE95 ~ URBLEVEL, data = ap)
Residuals:
Min 1Q Median 3Q Max
-14.873 -4.664 -1.174 3.639 37.569
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 6.3208 0.6251 10.11 <2e-16 ***
URBLEVEL 16.9877 1.0667 15.93 <2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 6.845 on 752 degrees of freedom
Multiple R-squared: 0.2522, Adjusted R-squared: 0.2512
F-statistic: 253.6 on 1 and 752 DF, p-value: < 2.2e-16
Regressao linear espacial - SAR
ap <- as.data.frame(cbind(crime_mg$MUNIC, crime_mg$AREA, crime_mg$INDICE94,
crime_mg$INDICE95, crime_mg$GINI_91, crime_mg$POP_94,
crime_mg$POP_RUR, crime_mg$POP_URB, crime_mg$POP_FEM,
crime_mg$POP_MAS, crime_mg$POP_TOT, crime_mg$URBLEVEL,
crime_mg$PIB_PC))
colnames(ap) <- c("ID",names(crime_mg@data[4:15]))
crime_mg_nb = poly2nb(crime_mg, queen=TRUE, row.names=crime_mg$X_COORD)
crime_mg_w <- nb2listw(crime_mg_nb, style="W")
sar.ap <- lagsarlm(crime_mg$INDICE95 ~ URBLEVEL,data=ap,crime_mg_w,method="eigen")
SARk_SSE <- sar.ap$SSE
SST <- sum((ap$INDICE95 - mean(ap$INDICE95))^2)
r2_SARk <- 1 - (SARk_SSE/SST)
r2_SARk
[1] 0.3314096
summary(sar.ap)
Call:lagsarlm(formula = crime_mg$INDICE95 ~ URBLEVEL, data = ap, listw = crime_mg_w, method = "eigen")
Residuals:
Min 1Q Median 3Q Max
-15.2482 -4.2371 -1.0771 3.3952 33.9250
Type: lag
Coefficients: (asymptotic standard errors)
Estimate Std. Error z value Pr(>|z|)
(Intercept) 2.35307 0.79725 2.9515 0.003163
URBLEVEL 13.93704 1.04857 13.2914 < 2.2e-16
Rho: 0.35437, LR test value: 66.163, p-value: 4.4409e-16
Asymptotic standard error: 0.044457
z-value: 7.9712, p-value: 1.5543e-15
Wald statistic: 63.539, p-value: 1.5543e-15
Log likelihood: -2486.1 for lag model
ML residual variance (sigma squared): 41.776, (sigma: 6.4634)
Number of observations: 754
Number of parameters estimated: 4
AIC: 4980.2, (AIC for lm: 5044.4)
LM test for residual autocorrelation
test value: 10.579, p-value: 0.001144
cat("Rˆ2 SAR: ",r2_SARk)
Rˆ2 SAR: 0.3314096
cat("Rˆ2 LM:",summary(lmK)$adj.r.squared)
Rˆ2 LM: 0.2511925
O modelo espacial SAR apresentou ganho de 8% versus o modelo linear simples.
Regressão linear simples:
ap <- as.data.frame(cbind(crime_mg$MUNIC, crime_mg$AREA, crime_mg$INDICE94,
crime_mg$INDICE95, crime_mg$GINI_91, crime_mg$POP_94,
crime_mg$POP_RUR, crime_mg$POP_URB, crime_mg$POP_FEM,
crime_mg$POP_MAS, crime_mg$POP_TOT, crime_mg$URBLEVEL,
crime_mg$PIB_PC))
colnames(ap) <- c("ID",names(crime_mg@data[4:15]))
head(ap)
lmK <- lm(formula = crime_mg$INDICE95 ~ URBLEVEL, data = ap)
summary(lmK)
Call:
lm(formula = crime_mg$INDICE95 ~ URBLEVEL, data = ap)
Residuals:
Min 1Q Median 3Q Max
-14.873 -4.664 -1.174 3.639 37.569
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 6.3208 0.6251 10.11 <2e-16 ***
URBLEVEL 16.9877 1.0667 15.93 <2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 6.845 on 752 degrees of freedom
Multiple R-squared: 0.2522, Adjusted R-squared: 0.2512
F-statistic: 253.6 on 1 and 752 DF, p-value: < 2.2e-16
Regressão linear GWR:
coords <- cbind(crime_mg$X_COORD,crime_mg$Y_COORD)
colnames(coords) <- c("X","Y")
ap <- as.data.frame(cbind(crime_mg$MUNIC, crime_mg$AREA, crime_mg$INDICE94,
crime_mg$INDICE95, crime_mg$GINI_91, crime_mg$POP_94,
crime_mg$POP_RUR, crime_mg$POP_URB, crime_mg$POP_FEM,
crime_mg$POP_MAS, crime_mg$POP_TOT, crime_mg$URBLEVEL,
crime_mg$PIB_PC))
colnames(ap) <- c("ID",names(crime_mg@data[4:15]))
bwGauss <- gwr.sel(crime_mg$INDICE95 ~ URBLEVEL,data=ap,coords=coords,adapt=TRUE,method="aic",
gweight=gwr.Gauss,verbose=FALSE)
gwr.ap <- gwr(crime_mg$INDICE95 ~ URBLEVEL, data=ap,coords=coords,bandwidth=bwGauss,
gweight=gwr.Gauss,adapt=bwGauss,hatmatrix=TRUE)
gwr.ap
Call:
gwr(formula = crime_mg$INDICE95 ~ URBLEVEL, data = ap, coords = coords,
bandwidth = bwGauss, gweight = gwr.Gauss, adapt = bwGauss,
hatmatrix = TRUE)
Kernel function: gwr.Gauss
Adaptive quantile: 0.013237 (about 9 of 754 data points)
Summary of GWR coefficient estimates at data points:
Min. 1st Qu. Median 3rd Qu. Max. Global
X.Intercept. -4.5982 4.2098 6.4983 9.3427 29.7712 6.3208
URBLEVEL -8.3511 10.0699 15.4027 20.0811 39.7302 16.9877
Number of data points: 754
Effective number of parameters (residual: 2traceS - traceS'S): 102.7988
Effective degrees of freedom (residual: 2traceS - traceS'S): 651.2012
Sigma (residual: 2traceS - traceS'S): 6.127271
Effective number of parameters (model: traceS): 71.52246
Effective degrees of freedom (model: traceS): 682.4775
Sigma (model: traceS): 5.985225
Sigma (ML): 5.694282
AICc (GWR p. 61, eq 2.33; p. 96, eq. 4.21): 4923.585
AIC (GWR p. 96, eq. 4.22): 4834.391
Residual sum of squares: 24448.34
Quasi-global R2: 0.4810664
GWR_SSE <- gwr.ap$results$rss
r2_GWR <- 1 - (GWR_SSE/SST)
r2_GWR
[1] 0.4810664
cat("Coeficientes LM:",summary(lmK)$coefficients)
Coeficientes LM: 6.320829 16.98768 0.6250511 1.066745 10.1125 15.92479 1.258985e-22 2.029188e-49
cat("Rˆ2 LM:",summary(lmK)$adj.r.squared)
Rˆ2 LM: 0.2511925
cat("Coeficientes GWR:",summary(gwr.ap$lm$coefficients))
Coeficientes GWR: 6.320829 8.987542 11.65426 11.65426 14.32097 16.98768
cat("Rˆ2 GWR:",r2_GWR)
Rˆ2 GWR: 0.4810664
cat("Rˆ2 LM:",summary(lmK)$adj.r.squared)
Rˆ2 LM: 0.2511925
cat("Rˆ2 SAR: ",r2_SARk)
Rˆ2 SAR: 0.3314096
cat("Rˆ2 GWR: ",r2_GWR)
Rˆ2 GWR: 0.4810664
Sim, GWR Aumenta 15% o ganho em relação ao SAR, que já era 8% maior que o modelo simples.
Implementação do modelo multivariado stepwise - regressão simples:
ap <- as.data.frame(cbind(crime_mg$MUNIC, crime_mg$AREA, crime_mg$INDICE94,
crime_mg$INDICE95, crime_mg$GINI_91, crime_mg$POP_94,
crime_mg$POP_RUR, crime_mg$POP_URB, crime_mg$POP_FEM,
crime_mg$POP_MAS, crime_mg$POP_TOT, crime_mg$URBLEVEL,
crime_mg$PIB_PC))
colnames(ap) <- c("ID",names(crime_mg@data[4:15]))
lm.ap <- step(lm(crime_mg$INDICE95 ~ ., data=ap))
Start: AIC=2234.09
crime_mg$INDICE95 ~ ID + AREA + INDICE94 + GINI_91 + POP_94 +
POP_RUR + POP_URB + POP_FEM + POP_MAS + POP_TOT + URBLEVEL +
PIB_PC
Df Sum of Sq RSS AIC
- AREA 1 3.8 14103 2232.3
- ID 1 6.2 14106 2232.4
- POP_MAS 1 10.7 14110 2232.7
- POP_FEM 1 14.5 14114 2232.9
- POP_RUR 1 15.7 14115 2232.9
- POP_URB 1 18.3 14118 2233.1
- POP_TOT 1 18.3 14118 2233.1
- POP_94 1 18.6 14118 2233.1
<none> 14100 2234.1
- PIB_PC 1 41.8 14141 2234.3
- GINI_91 1 232.1 14332 2244.4
- URBLEVEL 1 480.7 14580 2257.4
- INDICE94 1 18251.9 32351 2858.3
Step: AIC=2232.29
crime_mg$INDICE95 ~ ID + INDICE94 + GINI_91 + POP_94 + POP_RUR +
POP_URB + POP_FEM + POP_MAS + POP_TOT + URBLEVEL + PIB_PC
Df Sum of Sq RSS AIC
- ID 1 6.4 14110 2230.6
- POP_MAS 1 11.7 14115 2230.9
- POP_FEM 1 15.3 14119 2231.1
- POP_RUR 1 16.6 14120 2231.2
- POP_URB 1 18.5 14122 2231.3
- POP_TOT 1 20.5 14124 2231.4
- POP_94 1 20.7 14124 2231.4
<none> 14103 2232.3
- PIB_PC 1 41.4 14145 2232.5
- GINI_91 1 231.5 14335 2242.6
- URBLEVEL 1 477.9 14581 2255.4
- INDICE94 1 18434.4 32538 2860.6
Step: AIC=2230.64
crime_mg$INDICE95 ~ INDICE94 + GINI_91 + POP_94 + POP_RUR + POP_URB +
POP_FEM + POP_MAS + POP_TOT + URBLEVEL + PIB_PC
Df Sum of Sq RSS AIC
- POP_MAS 1 12.7 14122 2229.3
- POP_RUR 1 15.9 14126 2229.5
- POP_FEM 1 16.2 14126 2229.5
- POP_URB 1 18.4 14128 2229.6
- POP_TOT 1 19.5 14129 2229.7
- POP_94 1 19.7 14129 2229.7
<none> 14110 2230.6
- PIB_PC 1 41.1 14151 2230.8
- GINI_91 1 228.9 14339 2240.8
- URBLEVEL 1 474.5 14584 2253.6
- INDICE94 1 18464.5 32574 2859.5
Step: AIC=2229.31
crime_mg$INDICE95 ~ INDICE94 + GINI_91 + POP_94 + POP_RUR + POP_URB +
POP_FEM + POP_TOT + URBLEVEL + PIB_PC
Df Sum of Sq RSS AIC
- POP_FEM 1 3.5 14126 2227.5
- POP_RUR 1 13.9 14136 2228.1
- POP_TOT 1 19.2 14142 2228.3
- POP_URB 1 19.4 14142 2228.3
- POP_94 1 19.6 14142 2228.4
<none> 14122 2229.3
- PIB_PC 1 40.7 14163 2229.5
- GINI_91 1 230.6 14353 2239.5
- URBLEVEL 1 474.7 14597 2252.2
- INDICE94 1 18477.9 32600 2858.1
Step: AIC=2227.5
crime_mg$INDICE95 ~ INDICE94 + GINI_91 + POP_94 + POP_RUR + POP_URB +
POP_TOT + URBLEVEL + PIB_PC
Df Sum of Sq RSS AIC
- POP_URB 1 16.2 14142 2226.4
- POP_TOT 1 18.9 14145 2226.5
- POP_94 1 19.3 14145 2226.5
- POP_RUR 1 21.3 14147 2226.6
<none> 14126 2227.5
- PIB_PC 1 39.0 14165 2227.6
- GINI_91 1 227.1 14353 2237.5
- URBLEVEL 1 476.7 14603 2250.5
- INDICE94 1 18547.7 32674 2857.8
Step: AIC=2226.37
crime_mg$INDICE95 ~ INDICE94 + GINI_91 + POP_94 + POP_RUR + POP_TOT +
URBLEVEL + PIB_PC
Df Sum of Sq RSS AIC
- POP_TOT 1 18.2 14160 2225.3
- POP_94 1 18.7 14161 2225.4
- POP_RUR 1 21.7 14164 2225.5
<none> 14142 2226.4
- PIB_PC 1 39.4 14181 2226.5
- GINI_91 1 263.6 14406 2238.3
- URBLEVEL 1 465.1 14607 2248.8
- INDICE94 1 18809.7 32952 2862.2
Step: AIC=2225.34
crime_mg$INDICE95 ~ INDICE94 + GINI_91 + POP_94 + POP_RUR + URBLEVEL +
PIB_PC
Df Sum of Sq RSS AIC
- POP_94 1 2.1 14162 2223.4
- POP_RUR 1 21.6 14182 2224.5
- PIB_PC 1 37.4 14198 2225.3
<none> 14160 2225.3
- GINI_91 1 257.1 14417 2236.9
- URBLEVEL 1 448.5 14609 2246.8
- INDICE94 1 18805.2 32966 2860.5
Step: AIC=2223.45
crime_mg$INDICE95 ~ INDICE94 + GINI_91 + POP_RUR + URBLEVEL +
PIB_PC
Df Sum of Sq RSS AIC
- POP_RUR 1 21.0 14183 2222.6
- PIB_PC 1 36.1 14199 2223.4
<none> 14162 2223.4
- GINI_91 1 255.2 14418 2234.9
- URBLEVEL 1 447.1 14610 2244.9
- INDICE94 1 19048.5 33211 2864.1
Step: AIC=2222.57
crime_mg$INDICE95 ~ INDICE94 + GINI_91 + URBLEVEL + PIB_PC
Df Sum of Sq RSS AIC
- PIB_PC 1 36.5 14220 2222.5
<none> 14183 2222.6
- GINI_91 1 234.5 14418 2232.9
- URBLEVEL 1 454.0 14637 2244.3
- INDICE94 1 19027.5 33211 2862.1
Step: AIC=2222.5
crime_mg$INDICE95 ~ INDICE94 + GINI_91 + URBLEVEL
Df Sum of Sq RSS AIC
<none> 14220 2222.5
- GINI_91 1 252.2 14472 2233.8
- URBLEVEL 1 550.7 14771 2249.2
- INDICE94 1 19206.1 33426 2864.9
lm.ap
Call:
lm(formula = crime_mg$INDICE95 ~ INDICE94 + GINI_91 + URBLEVEL,
data = ap)
Coefficients:
(Intercept) INDICE94 GINI_91 URBLEVEL
4.6200 0.8295 -5.8052 5.3346
summary(lm.ap)
Call:
lm(formula = crime_mg$INDICE95 ~ INDICE94 + GINI_91 + URBLEVEL,
data = ap)
Residuals:
Min 1Q Median 3Q Max
-16.6372 -2.4362 -0.2178 2.3819 29.1584
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 4.62004 0.72718 6.353 3.65e-10 ***
INDICE94 0.82950 0.02606 31.827 < 2e-16 ***
GINI_91 -5.80519 1.59168 -3.647 0.000283 ***
URBLEVEL 5.33462 0.98980 5.390 9.47e-08 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 4.354 on 750 degrees of freedom
Multiple R-squared: 0.6982, Adjusted R-squared: 0.697
F-statistic: 578.3 on 3 and 750 DF, p-value: < 2.2e-16
lmKmv <- lm(formula = crime_mg$INDICE95 ~ INDICE94 + GINI_91 + URBLEVEL,
data = ap)
summary(lmKmv)
Call:
lm(formula = crime_mg$INDICE95 ~ INDICE94 + GINI_91 + URBLEVEL,
data = ap)
Residuals:
Min 1Q Median 3Q Max
-16.6372 -2.4362 -0.2178 2.3819 29.1584
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 4.62004 0.72718 6.353 3.65e-10 ***
INDICE94 0.82950 0.02606 31.827 < 2e-16 ***
GINI_91 -5.80519 1.59168 -3.647 0.000283 ***
URBLEVEL 5.33462 0.98980 5.390 9.47e-08 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 4.354 on 750 degrees of freedom
Multiple R-squared: 0.6982, Adjusted R-squared: 0.697
F-statistic: 578.3 on 3 and 750 DF, p-value: < 2.2e-16
Implementação do modelo multivariado stepwise - regressão SAR
ap <- as.data.frame(cbind(crime_mg$MUNIC, crime_mg$AREA, crime_mg$INDICE94,
crime_mg$INDICE95, crime_mg$GINI_91, crime_mg$POP_94,
crime_mg$POP_RUR, crime_mg$POP_URB, crime_mg$POP_FEM,
crime_mg$POP_MAS, crime_mg$POP_TOT, crime_mg$URBLEVEL,
crime_mg$PIB_PC))
colnames(ap) <- c("ID",names(crime_mg@data[4:15]))
sarmv.ap <- lagsarlm(crime_mg$INDICE95 ~ INDICE94 + GINI_91 +
URBLEVEL,data=ap,crime_mg_w,method="eigen")
SARkmv_SSE <- sarmv.ap$SSE
SST <- sum((ap$INDICE95 - mean(ap$INDICE95))^2)
r2_SARkmv <- 1 - (SARkmv_SSE/SST)
r2_SARkmv
[1] 0.7029335
summary(sarmv.ap)
Call:lagsarlm(formula = crime_mg$INDICE95 ~ INDICE94 + GINI_91 + URBLEVEL,
data = ap, listw = crime_mg_w, method = "eigen")
Residuals:
Min 1Q Median 3Q Max
-16.43661 -2.38664 -0.19365 2.28926 28.43605
Type: lag
Coefficients: (asymptotic standard errors)
Estimate Std. Error z value Pr(>|z|)
(Intercept) 3.36490 0.82039 4.1016 4.103e-05
INDICE94 0.80411 0.02706 29.7163 < 2.2e-16
GINI_91 -5.30895 1.58190 -3.3561 0.0007906
URBLEVEL 4.67297 0.99603 4.6916 2.711e-06
Rho: 0.10647, LR test value: 10.464, p-value: 0.0012169
Asymptotic standard error: 0.033258
z-value: 3.2013, p-value: 0.0013679
Wald statistic: 10.249, p-value: 0.0013679
Log likelihood: -2171.899 for lag model
ML residual variance (sigma squared): 18.562, (sigma: 4.3083)
Number of observations: 754
Number of parameters estimated: 6
AIC: 4355.8, (AIC for lm: 4364.3)
LM test for residual autocorrelation
test value: 0.10236, p-value: 0.74902
Comparação dos modelos Multi-variados de SAR e LM:
cat("Rˆ2 SAR: ",r2_SARk, "\n")
Rˆ2 SAR: 0.3314096
cat("Rˆ2 LM:",summary(lmK)$adj.r.squared,"\n")
Rˆ2 LM: 0.2511925
cat("Rˆ2 SAR MV: ",r2_SARkmv, "\n")
Rˆ2 SAR MV: 0.7029335
cat("Rˆ2 LM MV:",summary(lmKmv)$adj.r.squared,"\n")
Rˆ2 LM MV: 0.6969645
No modelo multivariado o ganho é de 1% sobre o modelo linear simples para o modelo SAR, porém o ganho foi de +37% com relação ao modelo original univariado.
Implementação do modelo multivariado stepwise - regressão GWR
bwGaussMV <- gwr.sel(crime_mg$INDICE95 ~ INDICE94 + GINI_91 + URBLEVEL,data=ap,coords=coords,adapt=TRUE,method="aic",
gweight=gwr.Gauss,verbose=FALSE)
gwrMV.ap <- gwr(crime_mg$INDICE95 ~ INDICE94 + GINI_91 +
URBLEVEL, data=ap,coords=coords,bandwidth=bwGauss,
gweight=gwr.Gauss,adapt=bwGaussMV,hatmatrix=TRUE)
gwrMV.ap
Call:
gwr(formula = crime_mg$INDICE95 ~ INDICE94 + GINI_91 + URBLEVEL,
data = ap, coords = coords, bandwidth = bwGauss, gweight = gwr.Gauss,
adapt = bwGaussMV, hatmatrix = TRUE)
Kernel function: gwr.Gauss
Adaptive quantile: 0.0411007 (about 30 of 754 data points)
Summary of GWR coefficient estimates at data points:
Min. 1st Qu. Median 3rd Qu. Max. Global
X.Intercept. -2.68291 2.58872 3.82700 8.52713 20.59749 4.6200
INDICE94 0.62641 0.75768 0.80116 0.83534 1.02374 0.8295
GINI_91 -27.22302 -11.11227 -5.72057 -0.53447 13.64553 -5.8052
URBLEVEL -2.24401 2.25196 5.03699 8.09010 12.19130 5.3346
Number of data points: 754
Effective number of parameters (residual: 2traceS - traceS'S): 66.21814
Effective degrees of freedom (residual: 2traceS - traceS'S): 687.7819
Sigma (residual: 2traceS - traceS'S): 4.135734
Effective number of parameters (model: traceS): 46.62188
Effective degrees of freedom (model: traceS): 707.3781
Sigma (model: traceS): 4.078046
Sigma (ML): 3.949956
AICc (GWR p. 61, eq 2.33; p. 96, eq. 4.21): 4313.115
AIC (GWR p. 96, eq. 4.22): 4257.927
Residual sum of squares: 11764.02
Quasi-global R2: 0.7503001
SST <- sum((ap$INDICE95 - mean(ap$INDICE95))^2)
GWR_MV_SSE <- gwrMV.ap$results$rss
r2_GWR_MV <- 1 - (GWR_MV_SSE/SST)
r2_GWR_MV
[1] 0.7503001
Comparando os resultados finais:
cat("Rˆ2 GWR: ",r2_GWR, "\n")
Rˆ2 GWR: 0.4810664
cat("Rˆ2 SAR: ",r2_SARk, "\n")
Rˆ2 SAR: 0.3314096
cat("Rˆ2 LM:",summary(lmK)$adj.r.squared,"\n")
Rˆ2 LM: 0.2511925
cat("Rˆ2 GWR MV: ",r2_GWR_MV, "\n")
Rˆ2 GWR MV: 0.7503001
cat("Rˆ2 SAR MV: ",r2_SARkmv, "\n")
Rˆ2 SAR MV: 0.7029335
cat("Rˆ2 LM MV:",summary(lmKmv)$adj.r.squared,"\n")
Rˆ2 LM MV: 0.6969645
Notamos que o modelo GWR tambem é beneficiado por uma análise multivariada, tendo aumentado 27%, passando de 48.1% para 75%.
ani_map <- tmap::tm_shape(crime_mg, simplify = 1) +
tmap::tm_fill() +
tmap::tm_shape(crime_mg) +
tmap::tm_fill(c("INDICE94","INDICE95")) +
tmap::tm_style(style = "natural", legend.outside = TRUE) +
tmap::tm_facets(free.scales.symbol.size = FALSE, nrow=1,ncol=1) +
tmap::tm_layout(main.title = "Evolucao do Indice de Criminalidade de 94-95") +
tmap::tm_polygons()
tmap_animation(ani_map, loop = TRUE, delay=200, filename = "CRIME_MG.gif")
Version: ImageMagick 7.0.8-10 Q16 x86_64 2018-08-15 https://www.imagemagick.org
Copyright: © 1999-2018 ImageMagick Studio LLC
License: https://www.imagemagick.org/script/license.php
Features: Cipher DPC HDRI Modules
Delegates (built-in): bzlib freetype jng jpeg ltdl lzma png tiff xml zlib
knitr::include_graphics("CRIME_MG.gif")